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Abstract 

Quasi-stationary distributions (QSD) have been widely studied since the pioneering work of 
Kolmogorov (1938), Yaglom (1947) and Sevastyanov (1951). They appear as a natural object when 
considering Markov processes that are certainly absorbed since they are invariant for the evolution 
of the distribution of the process conditioned on not being absorbed. They hence appropriately 
describe the state of the process at large times for non absorbed paths. Unlike invariant distributions 
for Markov processes, QSD are solutions of a non-linear equation and there can be 0, 1 or an infinity 
of them. Also, they cannot be obtained as Cesaro limits of Markovian dynamics. These facts make 
the computation of QSDs a nontrivial matter. We review different approximation methods for QSD 
that are useful for simulation purposes, mainly focused on Fleming- Viot dynamics. We also give 
some alternative proofs and extensions of known results. 
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1 Introduction 

Let Q = (q(x,y), x, y G A U {0}) be the rates matrix of a non-exploding continuous time 
Markov process Z(-,u), with initial condition u on a countable state space A U {0}. We 
assume that Q is irreducible in A and think of as an absorbing state. We also assume that 
absorption is certain. As usual we denote q(x, x) = — ^2 y ^ x q(x,y). 

Let be the absorption time defined by r M = inf{t > 0: Z(t,u) = 0} and T t \i the 
evolution of the process Z(-,u) conditioned on survival 

J2 z€ AKz)Pt(z,x) 



T t a{x) := P(Z(t,n) = x | r M > t) 
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Here Pt = e l ® is the semi-group associated with the transitions matrix Q. For fi such that 
fj,(z)\q(z, z)\ < oo, using the Kolmogorov equations for Z we obtain the Kolmogorov 
equations for the conditioned evolution 

yeA yeA (1.1) 

T (/i)(x) =fj,(x), x G A. 

In particular, we have the semi-group property T s+i = T s (T t ). We call (T 4 /i) the conditioned 
evolution. We can hence deduce from (11. ip that any limit point of T t \i is a fixed point of 
{T s , s > 0}. 

Definition 1.1. For a probability measure v in A, we say that it is a quasi-stationary 
distribution (QSD) if and only if for any t > we have T t v = v. 

Evaluating (II. ip at t — 0, we obtain a non-linear equation for v, i.e. v is a QSD if and 
only if it verifies 

° = J2 q ( y > x ) "(f) + °) "(f) "( x )- ( L2 ) 

Now, if v is a QSD, it follows from the Markov property that there exists 6 U > such that 

P(r u >t)= exp(-M)- 

Hence, a necessary condition for the existence of a QSD is the existence of a state x G A 
and # > such that E(e erX ) < oo. Here we slightly abuse notation and write t x instead 
of r Sx . Ferrari et. al (see below) showed in [27] that if A = N and for every t > 0, 
lmxj^oo P(t x < t) — 0, then this condition is also sufficient for the existence of a QSD (but 
does not imply uniqueness). Unfortunately, the additional condition lim^oo P(t x < t) = 
is not just technical (see [101 p. 9] for a simple counter-example). As a consequence, an 
absorbing process, even with integrable absorption times, might well have no QSD. It was 
also shown [301 SH 112] that there are a lot of cases where an infinite number of QSD do 
exist. These include for instance birth and death process with constant or linear rates and 
sub-critical branching processes. 

We now describe the dynamics of a process that is the cornerstone of several methods of 
both simulation and theoretical analysis of QSDs. 



The \i— return process. Given a probability /i in A, consider a Markov process with 
the following dynamics: it evolves according to the dynamics of Z until it hits 0, at which 
point, it instantaneously takes a new position in A with distribution fi. The transitions of 
this process are then given for x ^ y by: 



qt"(x, y) = q(x, y) + q(x, 0)/x(y). 



(1.3) 
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This process appears in the literature since the very beginning of the study of QSD 
([3 12]) but the approach was first discarded, arguing that there is a very weak relation 
between the QSD of Q and the invariant distribution of Z M and that this last one can be 
changed into almost any distribution by a suitable choice of \i. 

In [27], the authors considered the function that maps ji (-)■ where $(yu) is the 

invariant distribution of Z^. In that work, they realized that it was fruitful to study iterations 
of this map and they showed that v is a QSD for Q if and only if it is a fixed point of $. 
More importantly, they obtain the following theorem. 

Theorem 1.2 (Ferrari, Kesten, Martinez, Picco, [27]). Assume A = N and that P(t x > t) — >■ 
1 as x — > oo for every t > 0. Then a necessary and sufficient condition for the existence of 
a QSD is 

E(e 9rX ) < oo, for some x e N, 9 > 0. 

Their proof strategy is based on the fact that iterations of the map $ n (/i) are expected 
to converge to a QSD. 

Similar approaches were then exploited by different authors and with different goals. 

• The analysis of the dynamics of is the key elements for the proof of existence and 
uniqueness of QSD under Doeblin conditions given by Jacka and Roberts [37] . 

• Asselah and Castell [3] use the functional $ to prove existence of QSD for asymmetric 
conservative attractive particle systems in Z d . They work with Cesaro means of the 
iterations $ n (/x) rather than with the sequence itself. They prove that limit points of 
the Cesaro means are QSD and they also used this fact to prove properties (regularity) 
of the QSD relying on a priori bounds. 

• Aldous, Flannery and Palacios use also this idea to construct a method to simulate the 
QSD in finite spaces pQ. We give more details in Section [3j 

• Fleming- Viot dynamics are also linked to the dynamics of Z** . In fact, the distribution 
of a tagged particle of a stationary Fleming- Viot system with N particles follows the 
dynamics Z 11 , with fi being the empirical distribution of the process (in this case /i 
depends on time and is random and correlated with the tagged particle). When FV 
is in equilibrium, the tagged particle is expected to converge to Z v when A" goes to 
infinity, where v is a QSD. In Section [2TB1 we prove this fact for some models. 

• Also [21 E] deal with this map to approximate QSDs. 

As it can be seen, the /z-return process has motivated a lot of work in the field of QSDs 
in general and in particular in the problem of their simulation, which is non trivial since 
they can not be simulated directly, using for example the rejection method. In the sequel, 
we first review some results obtained in that direction in recent years. We mainly focus on 
Fleming- Viot (FV) dynamics, described in Section [21 as they constitute not only a method 
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to practically simulate the QSD and a theoretical tool to study QSDs, but also an interesting 
particle system raising many stimulating questions. 

We then provide a new proof for the convergence of Fleming- Viot dynamics to the con- 
ditioned evolution in compact time intervals and study the limiting behavior of a FV tagged 
particle. In Section [3j we explain the principles of two alternative methods for the compu- 
tation of the QSD. Section H] deals with related open problems. 

An updated bibliography on QSDs can be found at 

http : / /www .maths .uq.edu. au/ ~pkp /paper s/qsds/qsds . pdf 



2 Fleming- Viot dynamics 

We define the Fleming- Viot process driven by Z (or Q) as follows. It consists of N particles 
moving independently of the others as a continuous time Markov process with rates Q, until 
one attempts to jump to state 0. At that time, it immediately jumps back to A by choosing 
the position of one of the other particles chosen uniformly at random and starts afresh from 
this configuration. This process is Markov with state space A^ (here N is the number of 
particles) and its state at time t > is denoted £(£,•)• The variable £(t,i) denotes the 
position of the i-th particle at time t. The generator C N acts on functions / : A N -^las 
follows 

N ipN ^ 

= E E *) + 9(«0, 0) n^i^ ] (fi?'*) - /(£))> (2- 1 ) 

i=i xeA 

where ^ l ' x {i) = x, and for j ^ i, C' x (j) — £0)- We use {S^,t > 0} for the semi-group of 
FV, i.e. 

Sf/(£) :=£(/(£(*)) |£(0)=0 

The well definiteness of the process (given that the driving process is well defined) is 
immediate. Nevertheless we refer the reader to [UE1I29] or Section [2731 for (several) graphical 
constructions that are useful to infer different properties of the process. 

These dynamics were originally introduced by Burdzy, Ingemar, Holyst and March in the 
context of Brownian motion in a bounded domain in [18J, where it was conjectured that the 
empirical measure should converge to a deterministic solution of a PDE that represents the 
distribution of a Brownian particle conditioned on non-absorption. 

As it is described here, the process is a variation of the original measure-valued process 
introduced by Fleming and Viot [3T] (see also [2S])- Many alternative dynamics have been 
later proposed and studied. Their common feature is that all the particles evolve inde- 
pendently according to a common driving process until a selection mechanism takes place, 
replacing particles having specific positions or potentials. In Moran systems for example, 
particles branch according to a state-dependent rate and are replaced according to a fitness 
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functional of the positions. Brunet, Derrida and coauthors [T3J [T^l US1 US] on the other hand, 
introduced a system where particles evolve driven by one dimensional Brownian Motion and 
branch with a fixed rate, at which time the left-most particle is instantaneously eliminated. 
Although these systems are simpler to analyze than FV (in fact much finner results have 
been obtained for this model P, [39]), they can only be defined for random walks in Z + as 
driving process. 

Such models appear actually at the corner of several scientific fields. They were originally 
introduced in genetics to study the propagation of genes in a population. In a completely 
different context, they may represent the evolution of particles in a medium with an obstacle, 
and are linked with traveling waves equations in physics [13]. For more details, refer to [23] 
and the references therein. 

Coming back to FV, a generic question is hence to understand for a given driving Markov 
process, how the empirical density of the FV process at fixed times or under its stationary 
regime (if it exists) approximates the conditioned evolution or the QSD of the driving process 
respectively. This line of research was pursued by several authors and the conjecture was 
settled for many cases of interest where the driving process has a unique QSD, like for 
jump processes with finite state spaces [I], chains with Doeblin conditions [29] and diffusion 
processes on bounded domains [10J [191 E31 13"%] . 

An important and challenging problem then consists in studying the FV process when 
the driving process has an infinite number of QSDs. This question was recently answered 
for the case of a driving process given by sub-critical branching. In [5] it is shown that FV 
driven by this process is ergodic and that the empirical distribution of particles in equilibrium 
converges to the minimal QSD in the following sense. 

Definition 2.1. A QSD u* is called minimal if 



This phenomenon is called selection principle, in the sense that FV "selects" the minimal 
QSD among all the QSD. This fact also gives evidence in favor of the minimal QSD to be 
interpreted as the one that represents the (quasi) stationary behavior of the process before 
absorption takes place. 

Denote for x G A and £ G A N 



the proportion of ^-particles at position x. We also think of m(-, £) as a probability measure 
in A. With a slight abuse of notation, for a measure 7^ in A^ we write m(x, j N ) for the 
random measure m(x, £) where £ is a random configuration of particles with distribution 7 . 




M{E(t v ): v verifies (TP]) }. 




(2.2) 



i=l 
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Applying the generator of FV to the empirical distribution, we get 



C N m(y, •)(£) = ^2 l( x i y) m ( x i + 



N 



^2q(x,0)m(x,£)m(y,£,) 



(2.3) 



N-l 



which looks very similar to (11.11) . This provides another convincing heuristic for the approx- 
imation of the conditioned evolution and QSD by FV. 

It may be enlightening at this point to emphasize one of the main difficulties encountered 
when studying Fleming- Viot dynamics. Contrary to all the other mentioned genetic inter- 
acting jumps particle models, the FV process is not attractive: natural partial orders are not 
preserved by trajectories of the processes started from ordered initial positions. This compli- 
cates drastically the analysis and impedes the construction of simple coupling and stochastic 
comparisons for proving convergence properties. Another consistent difficulty arises because 
of the stochastic dependence between the times at which the selection mechanism takes place 
and the trajectories of the particles. As a consequence, the well definiteness of the process for 
all times when it is driven by diffusions is a difficult problenfl [TD1 El E3] • These difficulties 
can also be thought as somehow linked to the critical nature of those dynamics where the 
population size is kept constant. 

In the rest of this section, we aim at discussing: 

1. the convergence of the empirical distribution for fixed time intervals, 

2. the ergodicity and convergence of the empirical distribution under the stationary mea- 
sure of FV, 

3. the behavior of a tagged particle. 

We provide a new proof of convergence with order 1 / \/N for the first item and establish via 
a simple coupling argument the limiting distribution of a tagged particle. We also review 
the recent results obtained for long-time behaviors of the empirical distribution. 

2.1 Convergence to the conditioned evolution 

In this section, we investigate the convergence of the empirical distribution for fixed time 
intervals. In jl3], Villemonais proves a general result of convergence to the conditioned 
evolution with a rate of order 1 / \/N. The proof in [13] is based on martingale arguments 
under the assumption that the number of jumps of the FV process does not explode in 
compact time intervals. We here propose a new proof (leading to the same speed of conver- 
gence) working under quite general technical conditions on the transition rates and based 
on different techniques. It uses in particular the method of supersolutions. 

1 In countable spaces however, this follows immediately but this fact illustrates the possible complications 
that can arise with this dynamics. 
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Theorem 2.2. Assume J2 y eAl(y^ x ) — ^ f or a ^ x e A U {0}. If — J2 Z ^( z )q( z i z ) < °°> 
then 

sup E\m(x t t(t))-T t ii(z))\ < e 5CT + \\p - m(t(0))\\Tv) (2-4) 

o<t<T WJy J 



Remark 2.3. Note that the technical condition Y2 y eAl(y^ x ) < C is valid for a wide class of 
processes. In particular, we do not assume that rates are uniformly bounded nor the condition 

The proof of Theorem 12.21 is based on the following propositions. The first one is a 
useful decorrelation result proved in jl] and the second one proves the closeness between the 
semigroups associated with the conditional evolution of the driving process and the dynamics 
of the empirical measure of the FV process. 

Proposition 2.4 (Proposition 2 of jl]). Let Cq := sup xgA g(x, 0). For each t > 0, and any 
x, y G A 

O p 2C t 

sup \S?(m(x,-)m(y r )m-S t N m(x,-mS»m(y,-m]\ < — . (2.5) 



Remark 2.5. Proposition \2.J\ is proved in ^ in fact for processes with bounded rates, but 
the proof can be extended to this case with no difficulty. 

Proposition 2.6. Assume that Yl y eA ^{Vi x ) — ^ f or a ^ x e A U {0} . If — J2 Z ^( z )q( z j z ) < 
oo, then 

sup \Em(x,£{t)) - Ttn(x)\ < e 5CT ( 1 + - m(£(0))|| T v ] (2-6) 
Proposition 12.61 is based on the following lemma. 



Lemma 2.7. Assume that ]T^ga x ) < C for all x E AU {0}. Let e: [0, T] x A ->■ M be 

a function such that for x G A, t G [0, T] 

d 

— e(t, x) > ^2 q(z, x )e(t, z ) + c i(t, x)e(t, x) + c 2 (t, x) ^ q(z, 0)e(t, z), (2.7) 
zeA zeA 

and e(0,x) > for x G A, with c\ bounded on [0, T] and c 2 bounded and nonnegative. Then 
e(t,x) > for edit G [0,T], x G A. 

Proof. First, we can assume without loss of generality that C + c\ + c 2 < in [0, T]. If not, we 
consider w(t,x) := exp (— Xt)e(t, x), which verifies the same inequality with c\ replaced by 
c\ — A. Second, we can also assume that we have strict inequalities in (12. 7p . If not we consider 
e e (t, x) = e(t,x) + e which verifies (12.71) with strict inequalities thanks to C + c\ + c 2 < 0. 
Since e £ — > as e — > 0, if we prove e £ (t, x) > for all t, x, e, this implies our conclusion. 
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So, assuming that the conclusion of the lemma is false, there is a first time t* and a state 
x* such that e(t*,x*) = 0. At that time we have 

0> JU(f,<) > ^q(z,x*)e(t\y) + c 2 (f,x*)'%2q(z,0)e(t*,zf)>0, 

z^x* zeA 

a contradiction that concludes the proof. ■ 

Proof of Proposition 12.61 Define 

u(t, x) = St (m(x, •)) (£), and v(t, x) = T t u.(x) (2.8) 
and let <5(t, x) = u(t,x) — v(t,x). Recall that 
d 

—v(t, x) = j2 v( z i *m*> z ) + J2 M*> *m*' x )' ( 2 - 9 ) 

zeA zeA 

— u(i, x) = ^ q(z, x)u(t, z) + J2 0)u(i, *M*, ^) + a;, *), (2.10) 

zgA zeA 

where, 

= £ ?(*, 0) (^Y^f (Hz, -)m(x, •)) (0 - S^m(l, ■)(f)Sfm( ! c 1 -)(o) ■ 

z6A ^ ' 

(2.11) 

Proposition 12.41 implies that 

sup|i2(£;s,t)|<^^. (2.12) 

Then, for we get that: 

zeA zeA 

Call Cat := jr + — m(£(0))||Tv and consider the supersolution given by S(t, x) = 
5(t) := C N exp(5Ct). We have 

zGA zGA 

5(t) ^ x) + 5(t) £ g(z, 0)u(t, z) + 5(t)v(t, x) ^ 0) + R{£; x, t) < 

zeA zeA zeA 

2Ce 2C 0t ( 2 - 14 ) 

C N e 5Ct [ ?(*, x) + ?(*. 0)«(t, «) + s) £ q(z, 0) 1 



,zeA zeA zeA 



N 



< 



C N 5Ce 5Ct = ^-5(t, x) 
at 
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In the last line we have used that Cq < C and Cn > 1/N. Hence the function e(t,x) = 
S(t,x) —5(t,x) verifies (12. 7p with c\(t) = ^2 z& ^q(z,0)u(t, z) and C2(t) = v(t,x). Lemma [2~71 
then implies that 

6(t,x) < e 5Ct f 1 + - m(aO))\Wv) , x £ A, t e [0, T}. 

Proceeding in the same way with e(t, x) := 5(t, x) +5(t, x) instead of e we obtain the opposite 
inequality and hence 

\5(t,x)\ < e 5Ct (± + ||/i - m(£(0))|| TV J , x e A, t E [0,T]. (2.15) 



We can now proceed to the proof of the theorem. 

Proof of Theorem \2.2[ We apply the triangular inequality to get that for some k > and 
every t £ [0, T] 

E[\m(x,t(t))-T t u.{x)\] <E[\m(x,t(t))-u(t,x)\] + \u(t, x) - v(t, x)\. 
The control of the first term follows from Proposition 12.41 (with y = x) 

E[\m(x,C(t))-u(t,x)\] < E [|m(x,e(t)) -«(t,x)| 2 ] 1/2 

= \E{m{x^{t)?)-n{t,xf\ 1 / 2 
< ^/2/Ne Cot 

while the second term follows from (I2.15p . 



2.2 Long-time behavior 

The motivation for studying FV steams in particular from the exciting possibility of approx- 
imation of QSD, which are in general difficult objects to simulate. The kind of result one 
expects concerning the long-time behavior are 

1. For each N > 1, the Fleming- Viot process is ergodic, with invariant measure A . 

2. Let v* be the minimal QSD for the driving process, then 

lim f\m(x,£) - v*{x)\d\ N {£) = 0. (2.16) 
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A unique QSD. For Markov jump processes, these results have been proven in the case 
of finite state space in [I] and for processes satisfying a Doeblin condition in [29]. For finite 
state space, ergodicity is immediate while it follows from a "perfect simulation" argument in 
|29j . To prove the convergence of the empirical distribution under the invariant measure of 
FV is a much more involved issue. In both [U |29] , the key element of the proof is a control 
of the correlations between particles. In the case of a finite state space, the control of the 
correlations for a fixed time is obtained uniformly on the set of initial distributions (see !2.4p . 
This allows to extend this control under the stationary measure of FV. In [I], this step is 
done via a regeneration argument strongly connected to the Doeblin condition. 

An infinity of QSD: the Galton- Watson case. Many processes do actually have an 
infinity of quasi- stationary measures. The convergence mentioned in (I2.16P hence means 
that FV selects the minimal QSD. Proving these results becomes in this case much more 
challenging as one needs to prove not only tightness but also that the empirical measure of 
FV under the stationary measure belongs to the domain of attraction of the minimal QSD. 

An important example of processes owning an infinity of QSD is the subcritical branching 
process. For a full description of the number of QSD for branching processes see [H], and 
for general birth and death processes see [2D1 EDI EHl H2] • 

Up to our knowledge, the only result that deals with convergence of FV in equilibrium 
when the QSD is not unique is the following one, which is proved in [5]. 

Theorem 2.8 (Asselah, Ferrari, Groisman, Jonckheere, [5]). Consider a subcritical Galton- 
Watson process whose offspring law has some finite positive exponential moment. Let v* 
be the minimal quasi stationary distribution for the process conditioned on non- extinction. 



1. For each N > 1 the associated N -particle Fleming- Viot process is ergodic with invariant 



2. The empirical measure of FV converges to a deterministic measure which is the minimal 
QSD of the driving process: 



There are many other open questions concerning the long-time behavior of FV for large 
N. We deal with them in Section HI 

2.3 The tagged particle limit 

We now consider a particular particle (say particle 1) that we call the tagged particle. We 
show that under the invariant measure of FV, the tagged particle converges to a limiting 



Then 



measure A 




(2.17) 
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process that has the QSD distribution as marginal and is absorbed at exponential times. 
At those times, it jumps back instantaneously to A according to the QSD. A similar result 
holds out of equilibrium, but in this case the "return" jumps are according to the conditioned 
evolution. More precisely, the limiting process, (Y(t,{i),t > 0) is a continuous time non- 
homogeneous Markov jump process with initial condition \i and rates given by 

Qt(x,y) = q(x,y) + q(x,0)T t fi(y), x, y E A. (2.18) 

Observe that if /i is a QSD, then q^(x,y) is time-independent and also fi is invariant for 
(12.181) . Moreover, in this case (12.181) are the rates of the \x— return process described in the 
introduction. 

This problem has been previously addressed by Grigorescu and Kang [35] when the 
driving process is Brownian Motion in a bounded domain. 

In order to state our result we need to consider the Skorohod space D = D(IR + ,A) of 
right-continuous functions from M + — > A with finite left limits, equipped with the Skorohod 
topology [12]. We also introduce the notation £ 7 for a FV process with generator (12. ip and 
initial condition distributed according to 7. Finally, ^ is used when 7 = 5%. 

Proposition 2.9. Assume that Cq := sup a;eA q(x, 0) < 00 and let 7^ be a sequence of 
measures in A N such that 

lim sup E\m(x,C N (t)) -T t fi(x)\ = for all x £ A. (2.19) 

V-s-oo o<t<T 

Assume further that £ 7 (0, 1) converges in distribution to a random variable with distribu- 
tion \i. Then the tagged particle £ 7 (■, 1) converges in distribution to the process Y(-,/j,) as 
elements in D. 

As a consequence of Proposition 12.91 we obtain the following theorem. 

Theorem 2.10. Assume that for each N > 1, the FV process has an invariant measure X N 
and that for every x 

lim [ \m(x,£) - v(x)\d\ N (£) = 0. (2.20) 

AT->oo J 

Then, under X N the tagged particle converges in distribution to the stationary process Y(-, v). 

Remark 2.11. Hence Proposition ^. 9\ holds for every jump Markov process that verifies (I2.4p 
if m(x,£ 7 ) — > [i, which include a large class of absorbing processes as shown by Theorem 
\2.2l Theorem \2.10\ holds when one assumes in addition convergence of FV in equilibrium, 
which is a much more restricted class but includes the models of [29\/ . Of particular 
interest is the case of the driving process being a sub-critical branching process where an 
infinite number of QSD arise. In this case the results of /3]/ combined with Theorem \2.1U 
gives that the tagged particle in equilibrium converges to the stationary v* — return process, 
the probability v* being the minimal QSD. 
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2.3.1 Graphical construction 

The proof follows by coupling the trajectory of the limiting Markov process defined previously 
and the trajectory of a tagged particle of FV. To make the proof rigorous, we need to carefully 
build the coupling using a so called graphical construction. We first describe the construction 
for the FV process. 

Recall that C = sup xeA q(x,0) is the (maximum) absorption rate, and let 
g:=sup V" q{x,y); p{x,y) := q & y \ y -L x - p(x,x):=l- V] p(x,y). 

xGA Q 

xfcA yeA\{x} H yeA\{x) 

We assume for simplicity that q < oclE To each particle i, we associate two independent 
marked Poisson processes (oj(,u)Y) on R+, which we call respectively the internal and voter 
point processes, described as follows. 

• For each i, the internal Poisson process has intensity q. The marks are i.i.d uniform 
variables in [0, 1]. We use Uf for the mark of the k— th occurrence of uj. If the Poisson 
process rings at time t\ with a mark u\ whereas particle % is at x, then it jumps to y if 
and only if u\ G I(x, y), where for each x G A, {I(x, y) : y G A} is any partition of [0, 1] 
composed of Borel sets such that \I(x,y)\ = p(x,y). For details see [28j Chapter 1]. 

• The voter point process is also defined for each i. It has intensity Co and the marks 
are i.i.d. uniform variables on [0, 1] x [0, 1]. We denote them by (Bf, Vf). If uj rings 
at time t*, with a mark whereas particle % is at x, then it keeps its position if 
b\ > g ^°- > . If not, it jumps to zero and instantaneously comes back to A by jumping 
to y if and only if vf G Ji(t,y). Here {Ji(t,y) : y G A} is a (random) partition of [0, 1] 
by Borel sets such that \Ji(t,y) \ = ^ l{f (i, j) = y}. 

According to this procedure, particle i jumps to y with rate 

q — - — + Co Co iV _ 1 ^ mt,j) = Vh 

which coincides with the rates of (12.11) . 

We call oj = {{oj{ , toj), i G {1, . . . , N}) an i.i.d. sequence of marked PP. 

For any s < t real numbers, we denote by Ui[s,t) and Ui[s,t) the projections of the 
marked-times in the time period [s,t) and [s,t], respectively. 

We construct {£{t),t > 0} in such a way that £(t) is a function of the initial configuration 
^ and the time-marks a;[0, i], t > 0. Fix an initial configuration £ G A^, and t > 0. There 
is, almost surely, a finite number of time-marks within [0, £], say K\ let {s fc ,0 < k < K} be 
the ordered time-realizations with s° = 0. We build £(•) inductively as follows. 

2 this is not needed but otherwise once has to define a marked Poisson process for each potential position 
of a particle which complicates the exposition. 
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• At time s° = 0, the configuration is £(0) = £. 

• Assume £(s fc ) is known. For s £ [s k , s k+l ), set £(s) = £(s fc ). We describe now £(s fc+1 ). 

— If s k+1 corresponds to an internal time of particle i and mark u, we move particle 
% to y if and only if u £ J(£(s fc , i),y). 

— If corresponds to a voter time of particle i and mark (b,v), we move particle 
i to the position y if and only if b < 9 ^'°^ and u £ Jj(t, y). 

It is straightforward to check that {£(£),£ > 0}, as constructed above, has generator given 
by flUD, see [36]. 

2.3.2 The limiting process 

We now construct jointly with FV the limiting process {Y(t, //), t > 0} with initial condition 
/i as a deterministic function of (u[, u\) and a uniform random variable U. When is not 
necessary, we omit to display the dependence on \i and write just Y(-). Let {s k , < k < K} 
be the ordered time realizations of (u[,uY) in [0,t], with = 0. 

• In order to couple this process with FV, we use the same variable U uniformly dis- 
tributed in [0, 1] to define 

£(0,1) =j/ ^ Ue Ji(0,y) and F(0, //) = y ^ U £ J(0, y). 

The sets Ji(0, ?/) and J(0,y) are defined below. To get the correct distribution for 
the initial configuration of FV (7^) we construct the configuration of the remaining 
particles with the conditional distribution 7wU(o,i) (the procedure used to do this is not 
relevant). Next we describe the evolution of the limiting process. 

• Assume Y(s k ) is known. For s £ [s k , s k+1 ), set Y(s) = Y(s k ). We describe now Y(s k+1 ). 

— If s k+1 corresponds to an internal time of particle 1 and mark u, then Y jumps to 
y if and only if u £ I(Y(s k ),y). 

— If s k+l corresponds to a voter time of particle 1 and mark (b, v), then Y jumps to 
y if and only if b < and v £ J(t,y). Here {J(t,y): y £ A} is a partition of 
[0, 1] by Borel sets such that \J(t, y)\ = T t /i(y). 

It is also straight-forward to see that {Y(t),t > 0} has rates given by ( 12.181) . 

Observe that if v\ £ U weA [Ji(sf,j/) n J(s\,y)\ for every < k < K, then £(•, 1) = Y(-) 
in [0,t\. So we construct Jx(-,y) and J(-,y) in order to maximize their intersection. We 
consider a numbering of the elements of A (yi,l > 1) and define for t = 



^2^N(y r ),^2^N(y r 



r=l 



r=l 
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J(o,yi) 



1-1 



r=l 



r=l 



Here /i^v is the distribution of £(0, 1). For t > these sets take the form 

'i-i 



Ji{t,yi) 



\- J2i>l 1 {C( t , i ) =Vr} y> Ei>l i{f (*, *) = Vr} 



r=l 



i-1 



r=l 
/ 



iV - 1 



J(t,yi)= ^2 T tKyr),^2 T tKyr) 

_r=l r=l 

We also define 

^(t) := 1 - e U tf6A [J(*i,y) n J(fi,y)] for every fc with tj < *}. 

Note that ip(t) = 1 if the trajectories of both processes do not coincide on [0, £]. We have 
the following equation for the evolution of ip. 

Lemma 2.12. 



P(^(t) = l)<P(^(0) = l) + C o / E 



J2\rn(x,? N (s))-T s (ri(x)\ 



ds. 



(2.21) 



Proof. Since i/j(s—) is er([0, s))-measurable, we have 



and 



F(ifj(t) = 1) = P(V>(0) = 1) + / £ 



— P(^(s) = l|«7(a;[0,5))) 



ty>{s) = l\a(u>[0, 8 ))) 



ds, 



(2.22) 



j/GA 



(2.23) 



Observe that this also justifies the interchange of E and 4- in f)2.22p and completes the proof. 



Proof of Proposition \2. 9i . Recall that to prove convergence in D, it is enough to show con- 
vergence in D T = D([0, T],A), the Skorohod space of cadlag functions defined in [0, T] with 
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values in A, for every T > 0, cf. [121 Ch.3, Sec 16]. So, let Gt- D^ — )• M be a bounded, 
continuous function and define 

A N (T):= \Gr(t{;l))-G T (y{;n))\. 
Using that when ip(t) = 0, the trajectories of £(•, 1) and Y(-, fi) coincide on [0, t], we get that 
E(A N {T)) <E(A N {T)l{tff{T) = 0}) + E{A N (T)1{^{T) = 1}) 

^wGtW^p^t) = i) 

<2||G r || 00 [y;Mar)-^(a:)|+Cor sup V E(\m(x, ^(t)) - T tf i(x)\) ) . 

In the last line, we bound using Lemma 12.121 Now, since {T t fi, < t < T} is tight 
(as probability measures in A) and m(-,j ) — > fi, the convergence (I2.19P implies that 
{E(m(x,£ y (i))),0 < t < T, N > 1} is also tight and hence, for any e > 0, there exists a 
finite set A F such that 

sup ^n(x) + + E(m(x, ^ N '(£))) + E(T t u.(x)) < e. 



Next we use again the pointwise convergence (I2.19p for x £ A F and that /i^ — )• /i to get 



Kmswp E(A N (T)) < e. 



This concludes the proof. 



Proof of Theorem \2.1(A Since for each N, the measure X N is exchangeable and the empirical 
measure m(-,X N ) distributed according to X N converges to v in probability, we have that 
the tagged particle £(1) in equilibrium converges in distribution to v. So we can apply 
Proposition 12.91 with \i = v. ■ 



3 Other ways to approximate QSDs 

Since the QSD are eigenvalues of the transition matrices, in finite (small) spaces the compu- 
tation of the (unique) QSD can be addressed with standard (or not so standard, but generic) 
methods. However this can not be done (at least a priori) when the state space is huge or 
even countable. In this section we introduce two alternatives to FV dynamics to simulate 
the QSD. 
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3.1 The method of Aldous-Flannery-Palacios 

This method was proposed in [TJ to simulate the quasi-stationary distribution of a finite 
Markov chain in discrete time. To simplify the exposition, in this section we consider Z n as 
a discrete time Markov chain in a finite state space AU{0} with transition matrix P. We use 
P for the restriction of P to A and we assume P is irreducible and aperiodic. Hence, from 
Perron-Frobenius theorem we have the existence of a unique quasi- stationary probability 
measure v that verifies 

^2u(x)P{x,y) = Xu(y), ye A. (3.1) 

So, in order to compute v one has to solve (13. ip . but this is unpractical if A is large. Also 
the acceptance-rejection method is not useful since Pir > n) goes to zero exponentially fast. 

The method of Aldous-Flannery-Palacios is motivated on the /i— return process described 
in the introduction, i.e. the renewal process Z£ with transition matrix given by 

P»(x,y) = P(x,y) + P(x,0)n(y). 

That, is Z£ evolves as a Markov chain with transitions given by P but when it attempts to 
jump to zero, it comes back instantaneously to A according to \i. This is an ergodic Markov 
chain with invariant measure 3>(/f). If \i coincides with the QSD v then it is invariant for P v . 
Also the reciprocal holds, i.e.: v is a QSD if and only if Q{v) = v. This is the key fact used 
in [27] to prove their general existence result. The same fact is useful here to build a method 
approximating the QSD. The main idea is to replace v in the "return" process (which is only 
used to define the dynamics at absorption times) by the empirical distribution of the history 
of the process. 

Let M be the set of counting measures on A, and let (V n , fi n ) be the A x M— valued 
Markov chain with transitions defined by 

P(V n+1 = y, /i n+1 = fi + 8 v \V n = x, /i n = n) = P(x, y) + P(x, 0)^, (3.2) 

W 

where \fj,\ is the total mass of fi. In words, V n evolves as a Markov chain with transition matrix 
P, but when it attempts to jump to zero, it comes back to A instantaneously according to 
the empirical distribution of the history of the process fi. 

The following theorem was obtained by Aldous, Flannery and Palacios concerning this 
process. 

Theorem 3.1 (Aldous, Flannery and Palacios, [TJ). For the process defined by (13. 2 p we have 

lim — /i n = v. 

n— >oo fi 
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The proof follows by embedding fi n in the skeleton of a supercritical multitype branching 
process with offspring distribution given by L(x, y) = number of visits to the state y, for the 
process with transition matrix P, started at x. Details can be found in pp. 

This method, with improvements to gain efficiency was implemented by in [231 ESj to 
compute in particular the conditioned evolution and QSD of the contact process in different 
graphs. No proof of convergence are available in this case. 



3.2 Supercritical multitype branching 

The results of this section are the outcome of discussions with Pablo A. Ferrari. 

Consider the restricted matrix Q := (q(x, y) : x, y G A) and assume for the moment that 
A is finite. Let Q a be the matrix with entries 

q a (x, y) = q(x, y) + a5 xy , x, y G A. 

Here 8 xy is Kronecker's delta function. We consider a continuous-time multitype branching 
process as described by Athreya and Ney Chapter V.7]. The set of types is A. An 
individual of type x lives for a mean one exponential time and then splits into a random 
offspring distribution L(x) = (L(x,y), y G A) on N A with means given by m a (x,y) = 
q a (x, y) + S xy . We use X(t) = (X(t, x), x G A) for the number of individuals of each type at 
time t > 0. See [6] for details on this construction. It is well known (0, p. 202]) that the 
expected number of individuals of type y when we start the process with a unique individual 
of type x is given by E x (X(t, y)) = e tQa (x, y) for all x, y G A. 

Since A is finite, Perron- Frobenius theory (see also [22]) implies the existence of a principal 
eigenvalue X a for Q a (a real eigenvalue with real part larger than all other eigenvalues) and a 
positive left eigenvector u, normalized to be a probability measure in A. Observe that A a is 
an eigenvalue of Q a if and only if A Q = A + a, with A eigenvalue of Q and hence X a is positive 
if a is large enough. For those values of a, the associated multitype branching process with 
mean offspring matrix m Q is supercritical and the Kesten-Stigum theorem ([6], see also [32] 
for a conceptual proof) implies that 

hm - — — - = v, 



t— >oo 



almost surely on the event fl suvv := {|X(£)| > for all t > 0}. Now it is immediate that v 
is a left eigenvector for Q with eigenvalue A. Since ^ y( z\q{x,y) = —q(x,0), summing over 
y on both sides of 

^2v(x)q(x,y) = Xu(y), 

we get that A = —J2xeA h '( x ) c l( x ^) an d hence v verifies (11. 2p . which means that it is 
the unique QSD. This method is useful as a practical tool to simulate QSDs. It is valid 
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for finite A and in any other case where a Kesten-Stigum type theorem is available. So, 
for a given absorbing process with rates matrix Q, it is important to understand whether 
the supercritical multitype branching process described above behaves as predicted by the 
Kesten-Stigum theorem. 

4 Related open problems 

The different methods for simulating QSDs mentioned previously still raise many challenging 
open problems. We only mention here a few of them. 

Rates of convergence to equilibrium and error bounds. Very little information has 
been obtained concerning the invariant measure of FV and more generally about the rate of 
convergence and error bounds for the approximation of the QSD, for each of the mentioned 
methods. These would constitute very useful informations in particular for the FV dynamics: 
how fast is the FV process converging to its equilibrium for fixed N, and then how close is 
its empirical distribution from the QSD? 

The other eigenvectors. It has been shown that the FV dynamics allow the estimation of 
the maximal eigenvector of the infinitesimal generator associated with their driving process 
in many situations. One may naturally aims at showing similar results for other eigenvectors. 

In [18], the authors are interested in estimating eigenvectors of the Laplacian on bounded 
domains. They propose to adapt the dynamics of the FV process introducing different types 
of particles and considering as driving process random walks or Brownian motions. It would 
be of great interest to prove that this procedure converges in this case and can be adapted 
for other driving processes. 

Convergence of FV driven by drifted random walk. As it was pointed out before, 
the main questions related to FV are the ergodicity of the process for each N and the 
convergence of its empirical distribution. The case of FV driven by a drifted random walk is 
remarkably still an open problem. The difficulty stems from the fact that the driving process 
has an infinite number of QSDs together with the fact that it is not -R-positive (see |41j). 

The immortal particle for FV. Grigorescu and Kang [33] show that one can define the 
notion of immortal particle as the trajectory from which all particles proceed when time tends 
to infinity (in the context of diffusions processes). The convergence of the trajectory of the 
immortal particle towards a stationary regime (when it exists) has not been yet investigated. 
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Relation with N-BBM As alluded in the introduction, Brunet, Derrida and coauthors 
introduced a system where particles evolve driven by one dimensional Brownian Motion and 
branch with a fixed rate, at which time the left-most particle is instantaneously eliminated. 
Very precise quantitative informations have been investigated for these dynamics at various 
time scales [HI EH]- We conjecture that long time behaviors of this system properly renormal- 
ized are very similar to the one of the FV process. A systematic comparison of the empirical 
distributions and proof of convergence towards the minimal QSDs are open research avenues. 

Number of jumps of FV driven by diffusions in compact time intervals Though 
not included in the context of countable state space, we would like to mention the problem of 
the definition of the FV process for non-countable state space. As already underlined in the 
introduction, and defined as an open problem by K. Burdzy [17], it is not known in general 
under which conditions the number of jumps of the FV processes driven by diffusions on 
finite or infinite domain does not accumulate in compact times intervals. 

In [10] and [33], some conditions on the domain were shown to be sufficient to impede 
this phenomenon from happening. Also [11] partly deals with this problem. 
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